library(coefplot)
#install.packages("dotwhisker")
library(dotwhisker)
#install.packages("broom")
library(broom)
library(dplyr)
#install.packages("tidyr")
require(tidyr)
#install.packages("tidyverse")
require(tidyverse)
require(ggplot2)
require(stargazer)
#install.packages("huxtable")
require(huxtable)
library(Hmisc)
#install.packages("vctrs")
library(vctrs)
library(reshape2)
#install.packages("Rmisc")
library(Rmisc)
library(RColorBrewer)
#install.packages("wesanderson")
library(wesanderson)
#install.packages("doBy")
library(doBy)
#install.packages("binom")
library(binom)
#install.packages("miceadds")
library(miceadds)
library(gridExtra)
#install.packages("ggpubr")
library(ggpubr)
library(survey)
library(psych)
library(estimatr)
require(MASS)
require(Hmisc)
require(reshape2)
library(lmtest)
library("sandwich")
library(texreg)
library(interplot)
#devtools::install_github("strengejacke/strengejacke")
library(sjPlot)
#install.packages("ggeffects")
library(ggeffects)
#install.packages("ggcorrplot")
library("ggcorrplot")
#install.packages("Cairo")
library(Cairo)
#install.packages("scales")
library(scales)
#install.packages("lubridate")
library(lubridate)
require(multiwayvcov)
#install.packages("codebook")
#library(codebook)
#new_codebook_rmd()
#install.packages("future")
#library(future)
#install.packages("memisc")
#library(memisc)


setwd("~/II_Covid_Replication")

load(file = "bl.f.RData")


############
# Analyses #
############


### List of variables for different models, main analyses

mg.vars.b <- c("lost.job","cvd.hhs","cvd.cases.day")
mg.vars.e <- c(mg.vars.b,"cvd.gmb.vacc.rv","cvd.gmb.econ.rv","cvd.rch.vacc.rv","cvd.rch.econ.rv")
mg.vars.t <- c(mg.vars.e,"cvd.gmb.gvt")
mg.vars.m <- c(mg.vars.t,"mg.ab","rskavrs","nt.mg.out")
mg.vars.f <- c(mg.vars.m,"male","age","educ")

## Different baselines for job (Appendix C.1.)

full.nj <- mg.vars.f[mg.vars.f != "lost.job"]

mg.vars.bl0 <- c("lost.job",full.nj)
mg.vars.bl1 <- c("lost.job1",full.nj)
mg.vars.bl2 <- c("lost.job2",full.nj)

## Interaction between jobs and confidence in government (Main analyses)

mg.vars.int0 <- c("lost.job*cvd.gmb.gvt","lost.job",full.nj)


## List of all variable lists for OLS

mlist.ols <- list(mg.vars.b, mg.vars.e,mg.vars.t,mg.vars.m,mg.vars.f,mg.vars.bl0,mg.vars.bl1,mg.vars.bl2,mg.vars.int0)


## Object for storing regression results

sel.mg <- vector(length=length(mlist.ols), mode = "list")
names(sel.mg) <- c("b","beff","beffg","bgmig","bgmigdems","bgmigdemsj0","bgmigdemsj1",
                   "bgmigdemsj2","bgmigdemsint0")

## Running all OLS models

for (i in 1: length(mlist.ols))
{
  f <- as.formula(paste("mg.asp", 
                          paste(mlist.ols[[i]], collapse = " + "), 
                          sep = " ~ "))
  
  b.y.cl <- lm.cluster(f, cluster = "region.cln", data = bl.f) # Cluster standard errors by region

  
  b.y.x <- texreg::extract(b.y.cl) # For using with texreg
  sel.mg[[i]] <- b.y.x 
}


## Variables and labels for regression table

keep.vars <- list("lost.job1" = "Lost Job, CV-19\n(Baseline = No Job)", "lost.job2" = "Kept Job, CV-19\n(Baseline = No Job)", 
                  "lost.job10" = "No Job, CV-19\n(Baseline = Lost Job)", "lost.job12" = "Kept Job, CV-19\n(Baseline = Lost Job)", 
                  "lost.job20" = "No Job, CV-19\n(Baseline = Kept Job)", "lost.job21" = "Lost Job, CV-19\n(Baseline = Kept Job)", 
                  "cvd.hhs" = "HH Savings, CV-19", "cvd.gmb.vacc.rv" = "Effcy. CV-19 Vaccine, Gambia",
                  "cvd.gmb.econ.rv" = "Effcy. CV-19 Econ, Gambia", "cvd.rch.vacc.rv" = "Effcy. CV-19 Vaccine, Rich", 
                  "cvd.rch.econ.rv" = "Effcy. CV-19 Econ, Rich", "cvd.gmb.gvt" = "Trust Gov., CV-19",
                  "rskavrs" = "Risk Averse", "mg.ab" = "Ability to Migrate", 
                  "nt.mg.out" = "Migrant Networks", 
                  "male" = "Male", "age" = "Age", "educ" = "Education","cvd.cases.day" = "CV-19 Cases Cum. PMP",
                  "lost.job20:cvd.gmb.gvt" = "No Job, CV-19\n(Baseline = Kept Job):Trust Gov., CV-19",
                  "lost.job21:cvd.gmb.gvt" = "Lost Job, CV-19\n(Baseline = Kept Job):Trust Gov., CV-19",
                  "lost.job10:cvd.gmb.gvt" = "No Job, CV-19\n(Baseline = Lost Job):Trust Gov., CV-19",
                  "lost.job12:cvd.gmb.gvt" = "Kept Job, CV-19\n(Baseline = Lost Job):Trust Gov., CV-19",
                  "lost.job1:cvd.gmb.gvt" = "Lost Job, CV-19\n(Baseline = No Job):Trust Gov., CV-19",
                  "lost.job2:cvd.gmb.gvt" = "Kept Job, CV-19\n\n(Baseline = No Job):Trust Gov., CV-19")


###################################
############ Tables and Figures####
###################################


###########################
## Main Results (Table 2)#
###########################

texreg(list(sel.mg$b, sel.mg$beff, sel.mg$beffg,sel.mg$bgmigdems, sel.mg$bgmigdemsint0), custom.coef.map = keep.vars, 
       booktabs = T, dcolumn = T, stars = c(0.01, 0.05, 0.1), caption = "Main Analyses: Aspiration to Migrate", caption.above = T,
       custom.note = "SEs clustered by region")


#####################################
### Marginal Effects plot (Figure 3)#
#####################################

# running the model
int.blnojob.factor <- lm(mg.asp ~ factor(lost.job) + cvd.hhs + cvd.cases.day 
                         + cvd.gmb.vacc.rv + cvd.gmb.econ.rv + cvd.rch.vacc.rv + cvd.rch.econ.rv
                         + cvd.gmb.gvt + mg.ab + rskavrs + nt.mg.out
                         + male + age + educ + factor(lost.job)*cvd.gmb.gvt, data=bl.f)

##generating clustered standard errors
int.blnojob.factor.clust <- cluster.vcov(int.blnojob.factor, bl.f$region.cln)
int.blnojob.factor.mcse <- coeftest(int.blnojob.factor, int.blnojob.factor.clust)

int.blnojob.factor.mcse

## Pull out data for marginal effects plots
# coefficients
int.blnojob.factor.b2 <- int.blnojob.factor.mcse[2,1]
int.blnojob.factor.b3 <- int.blnojob.factor.mcse[10,1]
int.blnojob.factor.b22 <- int.blnojob.factor.mcse[17,1] # interaction term
# variances
int.blnojob.factor.varb2 <- int.blnojob.factor.clust[2,2]
int.blnojob.factor.varb3 <- int.blnojob.factor.clust[10,10]
int.blnojob.factor.varb22 <- int.blnojob.factor.clust[17,17]
# covariances
int.blnojob.factor.covb2b22 <- int.blnojob.factor.clust[2,17]
int.blnojob.factor.covb3b22 <- int.blnojob.factor.clust[10,17]

##transforming dataframe to graph    
int.blnojob.factor.df <- bl.f %>% mutate(cvd.gmb.gvt= seq(1, 5, length.out=829),
                                         int.blnojob.factor.conbx = int.blnojob.factor.b2+int.blnojob.factor.b22*cvd.gmb.gvt,
                                         int.blnojob.factor.consx = sqrt(int.blnojob.factor.varb2 + int.blnojob.factor.varb22*(cvd.gmb.gvt^2)+
                                                                           2*int.blnojob.factor.covb2b22*cvd.gmb.gvt),
                                         int.blnojob.factor.ax = 1.96*int.blnojob.factor.consx,
                                         int.blnojob.factor.upperx = int.blnojob.factor.conbx + int.blnojob.factor.ax,
                                         int.blnojob.factor.lowerx = int.blnojob.factor.conbx - int.blnojob.factor.ax)

# a rug

gmb.gvt.temp <- bl.f$cvd.gmb.gvt
a<-runif(829,min=0,max=0.6)
a<-a-0.2
gmb.gvt.jitter <-gmb.gvt.temp+a

blnojob.factor.plot <- ggplot(data = int.blnojob.factor.df) + 
  scale_y_continuous("Marginal Effect of Losing Job\n(Baseline: No Job)") + 
  scale_x_continuous(limits = c(1,5)) +
  geom_ribbon(aes(x=cvd.gmb.gvt, ymin=int.blnojob.factor.lowerx, ymax = int.blnojob.factor.upperx), alpha=.2) + 
  geom_line(aes(x=cvd.gmb.gvt, y=int.blnojob.factor.conbx, color="Lost Job"), size = 1)  +  
  geom_hline(yintercept = 0) +
  geom_rug(aes(x=gmb.gvt.jitter),sides="b") +
  scale_fill_manual(values = wes_palette(n=2, name="Darjeeling1"),
                    name = "",
                    labels = c("Lost Job")) + 
  xlab("Trust Gambian Government to Handle Pandemic") +
  ggtitle("Marginal Effect of Losing One's Job\n on the Aspiration to Migrate")+
  theme_bw(base_size = 12)+
  theme(text = element_text(size=20), 
        legend.position="none",
        plot.title = element_text(face="bold", hjust = 0.5)) 
blnojob.factor.plot  
ggsave("blnojob.factor.eps", width=20, height=12, unit="cm", dpi=300, device = cairo_ps)

####################################
# Different Baselines: Appendix C.1#
####################################

texreg(list(sel.mg$bgmigdemsj0,sel.mg$bgmigdemsj1,sel.mg$bgmigdemsj2), custom.coef.map = keep.vars, 
       booktabs = T, dcolumn = T, stars = c(0.01, 0.05, 0.1), caption = "Different Baselines: Aspiration to Migrate", caption.above = T)

##########################################
## Aspiration Ordinal Logit: Appendix C.2#
##########################################

f.ol <- as.formula(paste("mg.asp.f", paste(mg.vars.f, collapse = " + "), sep = " ~ "))
b.ol <- polr(f.ol, data = bl.f, Hess = T)
b.ol.cl <- coeftest(b.ol, vcov=vcovCL(b.ol, factor(bl.f$region.cln) ))


texreg(list(sel.mg$bgmigdems, b.ol.cl), custom.coef.map = keep.vars, 
       booktabs = T, dcolumn = T, stars = c(0.01, 0.05, 0.1), caption = "Ord Log: Aspiration to Migrate", caption.above = T,
       custom.model.names = c("OLS","Ord.L"))


#########################################
######### Descriptive Statistics ########
#########################################

###############################################
# Timeline Density Plots (Figure 1, main text)#
###############################################

cvd.att <- c("cvd.gmb.vacc", "cvd.rch.vacc","cvd.gmb.econ", "cvd.rch.econ")

cvd.df <- bl.f[,c("OLDUID",cvd.att)] 
long.cvd <- melt(cvd.df, id.vars = c("OLDUID"), measure.vars = cvd.att, value.name = "timelines", na.rm = T)

long.cvd.vacc <- subset(long.cvd, long.cvd$variable == "cvd.gmb.vacc" | long.cvd$variable == "cvd.rch.vacc" )
long.cvd.econ <- subset(long.cvd, long.cvd$variable == "cvd.gmb.econ" | long.cvd$variable == "cvd.rch.econ" )

pd <- position_dodge(0.2) # move them .05 to the left and right

p1 <- ggplot(data=long.cvd.vacc, aes(x=timelines, group=variable, color=variable, fill = variable)) +
  geom_density(adjust=2, alpha=0.02, na.rm = T, size = 0.7) + 
  ylab("") + labs(title="Vaccine Distribution") +
  scale_color_manual(values = wes_palette(n=2, name="Darjeeling1"),
                     name = "",
                     labels = c("The Gambia", "Rich Countries")) + 
  scale_fill_manual(values = wes_palette(n=2, name="Darjeeling1"),
                    name = "",
                    labels = c("The Gambia", "Rich Countries")) + 
  scale_x_discrete("", limits = c("Already", "Within year","Within 2 Years", "Distant future", "Never"), expand = c(0.03,0.03)) +
  scale_y_continuous(limits = c(0,1), expand = c(0,0)) +
  theme_bw(base_size = 17) + theme(plot.title = element_text(size = 17))


ggsave("Timeline_Vaccine.eps", width=22, height=10, unit="cm", device = cairo_ps)

p2 <- ggplot(data=long.cvd.econ, aes(x=timelines, group=variable, color=variable, fill = variable)) +
  geom_density(adjust=2, alpha=0.02, na.rm = T, size = 0.7) + 
  ylab("") + labs(title="Economic Recovery") +
  scale_color_manual(values = wes_palette(n=2, name="Darjeeling1"),
                     name = "",
                     labels = c("The Gambia", "Rich Countries")) + 
  scale_fill_manual(values = wes_palette(n=2, name="Darjeeling1"),
                    name = "",
                    labels = c("The Gambia", "Rich Countries")) + 
  scale_x_discrete("", limits = c("Already", "Within year","Within 2 Years", "Distant future", "Never"), expand = c(0.03,0.03)) +
  scale_y_continuous(limits = c(0,1.1), expand = c(0,0)) +
  theme_bw(base_size = 17) + theme(plot.title = element_text(size = 17))

ggsave("Timeline_Econ.eps", width=22, height=10, unit="cm", device = cairo_ps)


####################################
# Descriptives. Table 1, Main text #
####################################

# Proportions
lj <- table(bl.f$lost.job, useNA = "ifany")
male <- table(bl.f$male, useNA = "ifany")
reg <- table(bl.f$region.cln, useNA = "ifany")

# Others
describe(bl.f$cvd.hhs)
describe(bl.f$cvd.gmb.gvt)
describe(bl.f$rskavrs)
describe(bl.f$nt.mg.out)
describe(bl.f$age)
describe(bl.f$mg.asp)
describe(bl.f$mg.ab)
describe(bl.f$educ)
describe(bl.f$cvd.gmb.vacc.rv)
describe(bl.f$cvd.gmb.econ.rv)
describe(bl.f$cvd.rch.vacc.rv)
describe(bl.f$cvd.rch.econ.rv)


##################################
# Correlation matrix: Appendix C #
##################################

bl.covid.df <- bl.f[mg.vars.f]
dim(bl.f)
dim(bl.covid.df)
bl.covid.df$lost.job <- as.numeric(bl.covid.df$lost.job)


corr <- round(cor(bl.covid.df, use = "na.or.complete"), 2)
corr
colnames(corr) <- c("Lost Job", "HH Savings", "CV-19 Cases Cum", "Effcy. CV-19 Vaccine, Gambia",
                    "Effcy. CV-19 Econ, Gambia", "Effcy. CV-19 Vaccine, Rich", "Effcy. CV-19 Econ, Rich",
                    "Trust Gov., CV-19", "Ability to Migrate", "Risk Averse", "Migrant Networks", 
                    "Male", "Age", "Education")
rownames(corr) <- c("Lost Job", "HH Savings", "CV-19 Cases Cum", "Effcy. CV-19 Vaccine, Gambia",
                    "Effcy. CV-19 Econ, Gambia", "Effcy. CV-19 Vaccine, Rich", "Effcy. CV-19 Econ, Rich",
                    "Trust Gov., CV-19", "Ability to Migrate", "Risk Averse", "Migrant Networks", 
                    "Male", "Age", "Education")

ggcorrplot(corr,lab = TRUE, lab_size = 2, tl.cex = 7)

ggsave("corrmatrix.eps", width=12, height=10, unit="cm", device = cairo_ps)


########################################################
# Histograms: Figure 2 and Appendix E, Figures 3 and 4 #
########################################################

# Aspiration to migrate

aspiration.labs <- c("Dislike a\ngreat deal", "", "", "", "Like a\ngreat deal")

plot.aspiration <- ggplot(data = bl.f, aes(mg.asp)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)), width=0.25) +
  scale_y_continuous(labels = label_percent(accuracy=1L)) +
  scale_x_continuous(breaks=seq(1,5,1), labels = aspiration.labs) + 
  theme_bw(base_size = 12)+ 
  theme(text = element_text(size=20), 
        plot.title = element_text(face="bold", hjust = 0.5),
        axis.ticks.length = unit(.25, "cm"),
        axis.ticks.x = element_line(size=1)) +
  labs(x="", y="Percentage",
       title="Aspiration to Migrate",
       caption=str_wrap("Whether or not you think you would be ABLE to move abroad, how much would you LIKE to move to another country within the next year?", width=70)) 
plot.aspiration

ggsave("aspiration.eps", width=8, height=9, unit="in", device = cairo_pdf)

# Perceived ability to migrate

ability.labs <- c("Very Difficult", "", "", "", "Very Easy")

plot.ability <- ggplot(data = bl.f, aes(mg.ab)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)), 
           width=0.25) +
  scale_y_continuous(labels = label_percent(accuracy=1L)) +
  scale_x_continuous(breaks=seq(1,5,1), labels = ability.labs) + 
  theme_bw(base_size = 12)+ 
  theme(text = element_text(size=20), 
        plot.title = element_text(face="bold", hjust = 0.5),
        axis.ticks.length = unit(.25, "cm"),
        axis.ticks.x = element_line(size=1)) +
  labs(x="", y="Percentage",
       title="Perceived Ability to Migrate",
       caption=str_wrap("How easy do you think it would be for you, personally, to move abroad within the next year?", width=70)) 
plot.ability

ggsave("ability.eps", width=8, height=9, unit="in", device = cairo_pdf)

###########################################
# Histograms: Appendix E, Figures 5 and 6 #
###########################################

# Perceived ability to migrate through the backway

plot.ability.bw <- ggplot(data = bl.f, aes(mg.ab.bw)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)), width=0.25) +
  scale_y_continuous(labels = label_percent(accuracy=1L)) +
  scale_x_continuous(breaks=seq(1,5,1), labels = ability.labs) + 
  theme_bw(base_size = 12)+ 
  theme(text = element_text(size=20), 
        plot.title = element_text(face="bold", hjust = 0.5),
        axis.ticks.length = unit(.25, "cm"),
        axis.ticks.x = element_line(size=1)) +
  labs(x="", y="Percentage",
       title="Perceived Ability to Migrate\nThrough 'the Backway'",
       caption=str_wrap("Regardless of whether or not you would like to do this, how easy would it be for you to go through 'the backway' to Europe within the next year?", width=70)) 
plot.ability.bw

ggsave("ability_bw.eps", width=8, height=9, unit="in", device = cairo_pdf)

# Top destination
plot.dest <- ggplot(data = bl.f, aes(top.dest.num)) + 
  geom_bar(aes(y = (..count..)/sum(..count..)), width=0.25) +
  scale_y_continuous(labels = label_percent(accuracy=1L)) +
  scale_x_continuous(breaks=c(1,2),
                     labels = c("Within Africa", "Outside of Africa")) + 
  theme_bw(base_size = 12)+ 
  theme(text = element_text(size=20), 
        plot.title = element_text(face="bold", hjust = 0.5),
        axis.ticks.length = unit(.25, "cm"),
        axis.ticks.x = element_line(size=1)) +
  labs(x="", y="Percentage",
       title="Top Destination",
       caption=str_wrap("If you had to choose between your top destination within Africa and your top destination outside of Africa, which country do you think you'd be more likely to move to?", width=70)) 
plot.dest

ggsave("topdest.eps", width=8, height=9, unit="in", device = cairo_pdf)


#################################
### Covid Cases (Appendix C.3.)
################################

cvd.cases <- read.csv("owid-covid-data_gmb_jan2021.csv")
cvd.cases$date2 <- as.Date(parse_date_time(cvd.cases$date,"mdy"))

cvd1 <- cvd.cases[which(cvd.cases$date2 > as.Date("2021-01-06") & cvd.cases$date2 < as.Date("2021-01-28")),]

p<-ggplot(data=cvd1, aes(x=date2, y=new_cases)) +
  geom_bar(stat="identity") + labs(x= "Day", y="New Cases, Raw") + 
  scale_x_date(date_breaks = "1 day") + theme_bw(base_size = 12) +
  theme(axis.text.x = element_text(angle = 90)) 

q<-ggplot(data=cvd1, aes(x=date2, y=new_cases_per_million)) +
  geom_bar(stat="identity") + labs(x= "Day", y="New Cases, PMP") + 
  scale_x_date(date_breaks = "1 day") + theme_bw(base_size = 12) +
  theme(axis.text.x = element_text(angle = 90)) 

r<-ggplot(data=cvd1, aes(x=date2, y=total_cases)) +
  geom_bar(stat="identity") + labs(x= "Day", y="Cum. Cases, Raw") + 
  scale_x_date(date_breaks = "1 day") + theme_bw(base_size = 12) +
  theme(axis.text.x = element_text(angle = 90)) 


s<-ggplot(data=cvd1, aes(x=date2, y=total_cases_per_million)) +
  geom_bar(stat="identity") + labs(x= "Day", y="Cum. Cases, PMP") + 
  scale_x_date(date_breaks = "1 day") + theme_bw(base_size = 12) +
  theme(axis.text.x = element_text(angle = 90)) 


ggarrange(p, q, r,s, ncol=2, nrow=2)

ggsave("Cvd_Cases.pdf", width=8, height=9, unit="in", device = cairo_pdf)






